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Observations of the two gas giants show that both planets have dipolar magnetic fields: Jupiter's is very similar to Earth's 
magnetic field and Saturn's is very axisymmetric. Our main goal is to construct realistic numerical models that explain these 
features. 

While the small density jump across terrestrial iron cores allows to use the Boussinesq approximation, the picture is different 
for the gas giants. Here, the density decreases by a factor of around 5000 from the deep interior to the surface (1 bar level). 
Though most of this density jump is accommodated in the outer molecular envelopes, it may still be significant in the metallic 
dynamo region. Among other properties, the electrical conductivity also varies significantly with radius, being roughly constant 
in the metallic hydrogen region and decaying exponentially in the molecular envelope. We solve an anelastic numerical dynamo 
model (which differs from a fully compressible model by neglecting sound waves) to explore the effects of density stratification and 
electrical conductivity variation on magnetic field generation. 

We use an anelastic version of the MHD code MagIC with density jumps up to 245 and an electrical conductivity that decays 
exponentially in the outer 5-30% of the simulated shell. Previous simulations using constant conductivity showed that dipole- 
dominated magnetic fields are only found up to a density jump of 6. An increasing stratification progressively confines the most 
active convective region close to the outer boundary equator. Mean field models have shown that such a configuration prefers 
non-axisymmetric modes. The exponential conductivity decrease helps by separating magnetic field generation from the dominant 
convective region. 

For intermediate stratifications (6 < density jump < 148), the dipole component clearly dominates during short periods. Stable 
strongly dipolar solutions are found either when a large stratification (density jump > 148), at E — 10~ 4 , more clearly separates the 
dynamo from the dominant convective region, or when lower Ekman number are considered. 
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1. Introduction 

The dipolar magnetic fields of the gas giants, Jupiter and Sat- 
urn, are generated in their deeper interiors. Both planets mainly 
consist of a hydrogen-helium mixture. Due to the large pres- 
sures and temperatures reached inside these planets, hydrogen 
acquires metallic properties and becomes electrically conduct- 
ing (Chabrier et al., 1992; Fortney and Nettelmann, 2010). The 
transition happens at 85-90% of Jupiter's and 65% of Saturn's 
radii. The traditional view is that the lower metallic envelope 
likely hosts the dynamo of these planets, while the upper molec- 
ular envelope accommodates the observed fierce zonal jets. 
Higher densities, Lorentz forces and Ohmic diffusion would 
lead to a more sluggish dynamics in the metallic envelope and 
confine the zonal winds to the upper region. Traditional dynam- 
ical models therefore treat the two envelopes separately with 
dynamo simulations modeling only the metallic envelope and 
jet simulations concentrating on the molecular region. 
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The zonal jets have been investigated since the 70s by track- 
ing cloud features (see, for example, Ingersoll et al., 1979 for 
Jupiter and Sanchez Lavega, 1982 for Saturn). Their driving 
forces and depth are still debated. Some authors argue that 
they are a shallow weather phenomenon (Williams, 1978; Cho 
and Polvani, 1996) while others promote deeper-rooted jets that 
extend through the whole molecular envelope (Heimpel et al., 
2005; Jones and Kuzanyan, 2009; Gastine and Wicht, 2012). 
Gas giants emit roughly twice as much energy as they receive 
from the sun which argues for vigourous interior convection. 
In the rotationally-dominated dynamics ruling planetary atmo- 
spheres, interior convection quite naturally drives zonal winds 
via Reynolds stresses (i.e. a statistical correlation between the 
convective flow components; Christensen, 2002; Heimpel et al., 
2005). The winds follow a geostrophic structure, minimizing 
variations in the direction of the rotation axis, and therefore 
reach through the whole fluid atmosphere. Lian and Showman 
(2008) show that even when the forcing is restricted to a shallow 
weather layer the jets may reach much deeper into the planet. 
Kaspi et al. (2009), on the other hand, present an anelastic deep 
convection model where the equatorial jet is geostrophic and 
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the higher latitude jets are confined to the outer few percent in 
radius. 

Saturn's magnetic field is very axisymmetric and shows a 
strong concentration towards the poles (Cao et al., 2012) which 
seems incompatible with a classical Earth like dynamo model. 
A stably stratified layer at the top of the dynamo region (Chris- 
tensen and Wicht, 2008) or a completely different dynamo 
driven by differential rotation (Cao et al., 2012) are two alter- 
natives to explain these features. We will therefore concentrate 
on comparisons with Jupiter's field which is very similar to the 
geomagnetic field so that the well explored geodynamo mod- 
els also seem to apply at first sight. These typically adopt the 
Boussinesq approximation where the mild density stratification 
of Earth's core is simply ignored. Since this is certainly not 
justified in gas planets, some newer models therefore use the 
anelastic approximation which allows to incorporate the effects 
of the background density stratification while filtering out fast 
sound waves (Glatzmaier and Roberts, 1996; Stanley and Glatz- 
maier, 2010; Jones and Kuzanyan, 2009). In an extensive pa- 
rameter study, Gastine et al. (2012) show that dipole-dominated 
dynamos are rather rare in anelastic models once larger density 
stratifications are assumed. We further discuss this analysis in 
section 3. 

Ab initio calculations suggest that the electrical conductiv- 
ity of the molecular envelope can not be neglected altogether 
(Lorenzen et al., 2011; French et al., 2012). Due to the in- 
creasing degree of hydrogen ionization the conductivity expo- 
nentially increases with depth and matches the conductivity of 
the metallic region at the transition radius without any pro- 
nounced jump. Because of the high temperatures beyond the 
critical point no clear transitions is discernable. The classical 
separation of the dynamics for the two envelopes thus becomes 
questionable. Liu et al. (2008) argue that this has important 
consequences for the depth of the zonal winds. These should 
remain confined to a shallow outer layer where the conductiv- 
ity remains negligible. The strong shear associated with the 
zonal winds would otherwise create strong azimuthal magnetic 
field and lead to Ohmic heating incompatible with the observed 
luminosity (see however Glatzmaier, 2008). 

Stanley and Glatzmaier (2010) present an anelastic simula- 
tion of a relatively thin shell with exponentially decaying elec- 
trical conductivity to model the very outer part of the shell. 
Dipole-dominated magnetic field develops in the presence of 
strong zonal winds reaching through the simulated region. Ex- 
treme parameters (i.e. low Ekman and Prandtl number and 
high Rayleigh number) were used but a detailed discussion and 
systematic parameter study to disentangle the effects of den- 
sity stratification and varying conductivity are missing. Heim- 
pel and Gomez Perez (2011) also include a radial conductivity 
profile in their deep shell Boussinesq models, with a constant 
conductivity in the deeper interior and an exponential decay in 
the outer part. Once more well-pronounced deep-rooted zonal 
winds are compatible with dipole-dominated dynamo action. 

The present paper extends the work of Gastine et al. (2012), 
that we will refer to as GDW12 in the following, by adding 
an electrical conductivity profile allowing to incorporate both 
the metallic and the molecular envelope. Following Heimpel 



and Gomez Perez (2011), the profile is loosely based on the 
ab initio calculations from French et al. (2012) with a constant 
conductivity in the metallic and an exponentially decaying con- 
ductivity in the molecular region. The aim is to explore under 
which circumstances strong zonal surface winds and dipole- 
dominated dynamo action can coexist. 

In the next section, 2, we will describe the model with spe- 
cial attention to the anelastic formulation and the electrical con- 
ductivity profile. Section 3 presents the numerical results, first 
concentrating on the question of dipole-dominance and then on 
the dynamo mechanism. In the last section, 4, we summarize 
our main results and discuss their implication for the gas giants. 

2. Model 

2.1. Anelastic approximation 

The fluid and convective interior of the planet is modelled 
by solving the MHD equations in a rapidly-rotating spherical 
shell. In the past, the Boussinesq approximation was typically 
used, which neglects the background density and temperature 
variations. This is questionable in gas planets and, following 
Gilman and Glatzmaier (1981), Braginsky and Roberts (1995) 
and Lantz and Fan (1999), we therefore adopt the anelastic ap- 
proximation, where these background variations are now in- 
cluded but sound waves are ruled out by neglecting fast local 
density variations. 

We solve the equations in a dimensionless form (Christensen 
and Aubert, 2006; Schrinner et al., 2012), using the shell thick- 
ness as a lengthscale (d — r — r,, where r„ and r, are the 
outer and inner radii of the spherical shell, respectively) and 
the viscous diffusion time as a timescale (t v = d 2 /v, where v 
is the kinematic viscosity). Temperature and density are both 
non-dimensionalized by their values at the outer boundary, T„ 
and p . We employ constant entropy boundary conditions and 
use the imposed contrast across the shell as the entropy scale. 
There are no internal heat sources and all the heating coming 
into the shell via the inner boundary leaves it through the outer. 
While this is not the most realistic heating mode for gas gi- 
ants, it has been chosen to ease the comparisons with more 
classical Boussinesq simulations. The magnetic field is scaled 
by -Jp.AQ.po, where Aj is the magnetic diffusivity at r,-, p is the 
magnetic permeability and Q is the frequency of rotation of the 
shell. The radial magnetic diffusivity profile A(r) is normalized 
by its value at the inner boundary, A, = 1 /(cr,//), where cr ; is the 
electrical conductivity at r,-. The normalized magnetic diffusiv- 
ity is A = A(r)/Ai = l/(crp). 

The medium is assumed to be an electrically conducting 
ideal polytropic gas. Generally, dynamo simulations solve for 
small variations around an adiabatic hydrostatic background 
state that we mark with an overbar in the following. The tem- 
perature profile is then defined by the temperature gradient 
dT/dr = -g(r)/c p and the density profile by p(r) = T , where 
m is the polytropic index. We assume a simple gravity profile 
proportional to radius, which assumes density roughly homoge- 
neously distributed. The other extreme is to assume that all the 
mass is concentrated in the center, which leads to a gravity pro- 
file proportional to 1/r 2 (Gilman and Glatzmaier, 1981; Jones 
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et al., 2011). Gastine et al. (2012) have shown that both gravity 
profiles lead to very similar results. The true profile of the gas 
giants lies somewhere in-between. The temperature reference 
state is then given by 



7» = -c (^) +l+co, 
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N p = \n(pjlp ) is the number of density scale heights be- 
tween the inner and the outer boundaries of the shell and rj is 
the ratio between the corresponding radii (see Jones et al., 201 1 ; 
Gastine and Wicht, 2012, for the full derivation of the reference 
state). 

The dimensionless form of the anelastic equations is then 
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where the traceless rate-of-strain tensor with constant kinematic 
viscosity is given by 



S = 2p 



e, 7 - -SijV ■ u 



and 



1 / dui duj 

2 V dxj dxi 



), (8) 



where 6jj is the identity matrix. The viscous and ohmic heating 
contributions are given by 



Gv = 2p 
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The system of Eqs. (3-7) is governed by the dimensionless 
Ekman number, Rayleigh number, Prandtl number and mag- 
netic Prandtl number, respectively, 
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where the thermal diffusivity k, the specific heat c p and the 
value of gravity at the outer radius, g , are constant. 

2.2. Variable conductivity 

To simulate the variable electrical conductivity of hydrogen 
in the interior of Jupiter, we employ a profile that corresponds 
to an approximately constant conductivity in the metallic hy- 
drogen layer and an exponential decay in the outer molecular 
envelope. Both branches are matched via a polynomial that also 
ensures that the first radial derivative is continuous: 
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The exponential decay with a rate a starts at a radius r m where 
the normalized conductivity has decreased from cr,- = 1 to <x m . 
For convenience we also define the relative transition radius 
Xm - r m /r„ and we refer to it as a percentage. 

This kind of profile has first been used by Gomez-Perez et al. 
(2010) and it seems a fair first approximation to results from ab 
initio calculations (French et al., 2012). In our study, we mostly 
fix & m to 0.5 and a to 9 while \m is varied from 95% to 70%. 
Since the decay rate is much too small when compared with 
the ab initio results, we also include a few cases with a much 
steeper decay of a = 25. Despite the very slow conductivity 
decay at a — 9, our results show no drastic changes and thus 
prove that this value seems sufficient to capture the main effects. 

2.3. Numerical model 

For the numerical simulation of the model described above, 
we use the anelastic version of the MagIC code (Wicht, 2002; 
Gastine and Wicht, 2012). This is a pseudo-spectral code that 
solves Eqs. (3-7) in a spherical shell using a poloidal/toroidal 
decomposition of the vector fields pu and B, 
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where the scalar fields w, c are poloidal potentials and z, a 
are toroidal potentials. These potentials are then expanded in 
spherical harmonic functions up to degree and order l max in lat- 
itudinal and longitudinal directions (9 and cp, respectively) and 
in Chebyshev polynomials in the radial direction (up to degree 
N r ). 

For the parameter studies presented here, we use different 
resolutions, because the higher gradient in density and electri- 
cal conductivity demand an increase of both radial and hori- 
zontal resolution. The Chebyshev polynomials are truncated 



from N r = 73 to N r = 97 and l max ranges between 85 and 170. 
Each simulation ran for at least 1 magnetic diffusion time, ex- 
cept cases with lower Ekman number (see Tab. 2). 

The boundary conditions for the velocity of our models are 
no slip at the inner core boundary and stress-free at the outer 
boundary. The stress-free outer boundary condition is more 
appropriate for a gas giant and allows for a more dynamical 
behaviour when compared with a rigid boundary (Simitev and 
Busse, 2009; Schrinner et al., 2012; Gastine et al., 2012). The 
existence of an inner core in Jupiter and its possible size (if it 
exists) is still unclear. Here we assume a small and conducting 
solid inner core with rjr =0.2 (bigger than the model by Net- 
telmann et al., 2012, though), and we use a rigid flow boundary 
condition at the interface. 

Furthermore, we use constant entropy boundary conditions, 
which means that we don't consider internal heat sources in our 
models. We match the magnetic field to a diffusive solution at 
the inner boundary and to a potential field at the outer boundary. 

The Ekman number is either E = 10~ 4 or E — 10 5 for the sim- 
ulations presented here. The relative moderate value of E — 10 4 
allows to more extensively scan the other system parameters 
like Rayleigh number, density stratification N p and transition 
radius of the electrical conductivity profile;^,,,. At E- 10" 5 , we 
could only afford to run a few cases. We assume a Prandtl of 
Pr = 1 and magnetic Prandtl number at the inner boundary of 
Pm.i = 2 for E= 10~ 4 and Pm, = l for £= 10~ 5 . 
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Figure 1 : Dipolarity against the local Rossby number calculated only 
for the conducting region (Eq. 21). Colour corresponds to different 
transition radii: cyan —Xm = 70%, red - Xm = 80%, blue -Xm = 90%, 
green - Xm - 95%. The open symbols are cases with a no slip upper 
boundary, instead of stress-free. The error bars correspond to standard 
deviations of the time series of each case, for which the point itself is 
the time average. The cases in black color correspond to a = 25 and 
Xm = 80%. The three dashed lines connect three sets of cases for which 
we found two solutions, depending on the initial magnetic field. The 
cases indicated in dark yellow outlined in red (green) correspond to 
E = 10~ 5 simulations with^„, = 80% (x m = 95%). Finally the yellow 
symbol without the red line is the case from Heimpel and Gomez Perez 
(2011)of^ m = 80% and 77 = 0.35. 



Table 1 : Values of critical Rayleigh number (Ra cr ) and critical wave 
number (m cr ) for each N p . The values were obtained with a modified 
version of the linear code by Jones et al. (2009). 
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All together, we simulated 50 cases in a range of 4 to 46 
times supercritical. Gastine and Wicht (2012) examined anelas- 
tic convection for an aspect ratio of 0.6 and showed that the 
critical Rayleigh number increases with increasing stratification 
N p . Table 1 demonstrates that we observe a similar trend for the 
smaller aspect ratio 0.2 employed here. At a certain stratifica- 
tion the wave number jumps from lower to high values. This is 
the point where the onset of convection moves from the inner to 
the outer boundary (Gastine and Wicht, 2012). The transition 
is delayed for the lower Ekman number (see also Jones et al., 
2009). 

We consider a polytropic index of m-2 in our anelastic sim- 
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Figure 2: Dipolarity plotted against the ratio between the toroidal 
kinetic energy and the total kinetic energy, averaged in time and vol- 
ume. The symbols have the same definition as in Fig. 1. The three 
boxes mark the three different regimes discussed in the text. The three 
dashed lines connect three sets of cases for which we found two solu- 
tions, depending on the initial magnetic field. 



ulations. The explored density scale heights range from the 
Boussinesq case N p = to N p — 5.5 which corresponds to a 
density jump of p,/p - 245. Ab initio simulations suggest a 
Jovian stratification of N p = 8.5 (p,/p - 5000) from the bot- 
tom of the molecular hydrogen layer to the 1 bar level (Guillot, 
1999; French et al., 2012). However this is mostly concentrated 
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in the outermost 1% and our largest stratification indeed covers 
99% of the planet. 

2.4. Diagnostic parameters 

The parameters of all the numerical experiments discussed 
here are listed in Tab. 2 along with several quantities that char- 
acterize the solution. The amplitude of the zonal axisymmetric 
flow component is measured in terms of the Rossby number 

ROyfm . 
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and V is the volume of the spherical shell. The magnetic 
Reynolds number estimates the ratio of magnetic field produc- 
tion and diffusion. We use a modified form to account for the 
radial-dependent magnetic diffusivity, 
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The local Rossby number has been introduced by Chris- 
tensen and Aubert (2006) to quantify the relative importance 
of the advection term in the Navier-Stokes equation (Eq. 3): 
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Here, we use a modified version of Ro[ calculated only for the 
inner conducting region (r, < r < r m ), 
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with Uj the flow contribution of spherical harmonic degree I. 

The magnetic field strength is quantified by the Elsasser 
number which measures the ratio of Lorentz to Coriolis forces 
using the modified form 

A= — l - — f /— \r 2 dr. (23) 

J,., \p(r)A(r)l 

The geometry of the surface field is characterized by its dipo- 
larity 
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which measures the relative energy in the axial dipole contribu- 
tion at the outer boundary r a . Following Christensen and Aubert 
(2006), we restricted the magnetic field to spherical harmonic 
degrees and orders below 12 in Eq. 24. The time variability of 
diplarity is quantified by its standard deviation SDdip shown in 
Tab. 2. 

The properties described in this section are referred in Tab. 2 
by their time averages. 
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Figure 3: Time evolution of the dipolarity at the surface (solid black 
line) and the tilt angle of the dipole (dashed grey line) in degrees. The 
time is given in magnetic diffusion units. The parameters are: E = 10 -4 , 
N p = 3, Ra =4.3 Ra cr , X m = 80%. 
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Figure 4: Dipolarity against density stratification, for fixed transition 
radii: red - 0.8 r rj , blue - 0.9 r„, green - 0.95 r . The dashed lines simply 
gather the symbols with the same^,,,. The error bars are the same as in 
Fig. 1 . Thexm = 70% cases were excluded, since we didn't find dipolar 
solutions up to now. 



3. Results 

Before presenting the new simulations for a radial conductiv- 
ity profile, we briefly summarize the main findings of GDW12. 
They explored the effects of density stratification on numerical 
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N p = 0, X m = 95% N p = 4, X m = 95% 




Figure 5: Radial magnetic field at the outer boundary. The top row corresponds to^ m = 95% (cases 9d and 33 of Tab. 2, respectively) and the 
bottom row to Xm = 80% (cases 3 and 34 of the same table). The maps on the left are N p = cases and the maps on the right column are N p =4 
(top) and N p = 5 (bottom). Magnetic fields are given in units of Elsasser number. 



dynamo simulations with homogeneous electrical conductivity, 
focussing mainly on the question of magnetic field geometry. 
Ekman and Prandtl numbers are identical to most of the cases 
explored here. The stratification was varied up to N p — 3 and 
two different aspect ratios, rj = 0.2 and rj = 0.6, were explored 
with very similar results. 

For the smaller to intermediate stratifications N p < 2, only 
multipolar magnetic fields were found when the local Rossby 
number exceeds about Roi w 0. 1 This confirmed the findings of 
Christensen and Aubert (2006) who introduced the Rog to sep- 
arate dipole-dominated from multipolar dynamo simulations in 
their Boussinesq models. 

Below Roc m 0.1, GDW12 find a bistable regime where mul- 
tipolar and dipole-dominated solutions coexist and form two 
distinct branches. This bistability can likely be attributed to 
the fact that stress-free boundary conditions were employed 
(Simitev and Busse, 2009; Schrinner et al., 2012; Gastine et al., 
2012). For rigid conditions, only the dipolar branch seems to 
exist. Stress-free boundary conditions allow strong zonal winds 
to develop that compete with large scale magnetic fields. On the 
dipolar branch zonal winds are weak, on the multipolar branch 
they are stronger. 

For larger density stratifications, N p > 2, no dipole- 
dominated solutions were found. This was attributed to the fact 
that the focus of convective action moves progressively outward 
when the stratification is increased. When convective columns 
are mainly confined to a thin region close to the outer bound- 
ary, a new non-axisymmetric dynamo mode is preferred that 
has previously been observed in some mean field dynamo sim- 
ulations (Rudiger et al., 2003; Jiang and Wang, 2006). 



3.1. Dynamo regimes 

We start with discussing the different dynamo regimes based 
on the E = 10" 4 cases and come back to E = 10~ 5 solutions in 
section 3.3. Our new results for depth-dependent conductivity 
are subsumed in Figs. 1 and 2. Figure 1 displays the time aver- 
aged dipolarity at the outer boundary versus the time-averaged 
local Rossby number calculated for the metallic hydrogen re- 
gion only (r/r a < x,„)- 

Figure 2 shows the time-averaged dipolarity versus the rela- 
tive time-averaged zonal wind energy. Three distinct groups of 
solutions can be identified. 

When the weakly conducting layer is thin (Xm = 95% and 
Xm = 90%) we largely recover the results from GDW12: for 
weak stratification and Ro( < Ro[ c , there are two branches with 
either weak dipolarity and stronger zonal flows (cyan area in 
Fig. 2) or strong dipolarity and weaker zonal flows (green and 
blue symbols inside the yellow area in Fig. 2). Our results also 
agree with the findings by Christensen and Aubert (2006) in 
the sense that above a critical local Rossby number value Rof c 
only multipolar solutions can be found. However, the value of 
Ro c is significantly decreased from 0.1 to about 0.045. This is 
consistent with Gomez-Perez et al. (2010) who reported that in- 
troducing a thin weakly conducting layer in a dipole-dominated 
dynamo model resulted in a multipolar solution. The aspect ra- 
tio has an additional but milder effect on the transitional Rot 
value which is reduced for smaller inner cores. 

When the thickness of the weakly conducting layer is in- 
creased, we also find two distinct types of solutions: a new 
group of very time-dependent solutions (the pink area in Fig. 2) 
and a group of solutions with more stable strongly dipolar fields 
(within the yellow area). The latter results from combining 
a strong stratification (diamonds, N p = 5 and left-triangles, 
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Figure 6: Azimuthal averages of the zonal component of the flow. Each column of two plots has a different N p , namely 0, 1, 3, 4 and 5 from left 
to right. In the top row, the poloidal fieldlines are plotted on top of the zonal velocity contours, where solid lines are positive values and dashed 
are negative values. The dotted line corresponds to the transition r m = 80%. The bottom row shows the corresponding hydrodynamical solutions 
(i.e. without magnetic field). All cases have the same contour levels (the red corresponds to prograde and the blue to retrograde, the maximum 
and minimum values are 0.05 and -0.05, respectively, in units of Rossby number determined by Ro = u/(Slr„). 



N p = 5.5) with a thicker weak conductivity layer (red colored 
symbols, Xm - 80%). 

The new oscillatory solutions promoted by the more pro- 
nounced weak conductivity layer (red, Xm - 80% and cyan, 
Xm = 70%) are indicated by large error bars (within the pink 
box in Fig. 2), measuring the standard deviation around the 
time-averaged dipolarity. Figure 3 demonstrates that this strong 
time dependence results from the fact that these cases oscillate 
between dipolar and multipolar solutions. When the dipole is 
weak, it can significantly tilt away from the rotation axis and 
even undergo polarity reversals. 

3.2. Competition between zonal flows and magnetic field 

Figure 4 highlights the role of the density stratification at 
E = 10~ 4 . When the weakly conducting layer has a negligi- 
ble influence, dipole-dominated solutions can only be found for 
N p < 2, as well as bistability (Y shape of the green dashed 
line). For thicker and thus influential layers, however, the mean 
dipolarity increases with stratification. For low to intermediate 
N p , oscillatory cases, similar to the case displayed in Fig. 3, are 
found. For larger stratifications they finally settle into the very 
strongly dipolar and more stable solutions. 

Snapshots of the radial magnetic field at E = 10~ 4 at the 
outer boundary illustrate the different type of solutions in Fig. 5. 
The top row shows a dipole-dominated Boussinesq case and 



the multipolar solution caused by a strong density stratifica- 
tions (see GDW12). Note that the field presents a wavy struc- 
ture (m=l) typical for this multipolar dynamos with stress-free 
boundaries (Goudard and Dormy, 2008, GDW12). The weakly 
conducting layer is only thin here with Xm =95% and has negli- 
gible influence on the dynamo action. The second row depicts 
the new solution branches found for^„, = 80%. The left panel 
shows a snapshot of a Boussinesq multipolar case and the right 
panel illustrates the dipolar configuration found at strong strati- 
fications. At the equator, the dipole field is reduced in the N p -5 
cases. 

Why does the combination of a thicker weakly conduct- 
ing outer layer and a strong stratification help to get dipole- 
dominated solutions? We recall that GDW12 found two reasons 
why the dynamo process prefers to produce multipolar fields: 
for one, there is a competition between strong zonal jets and 
larger scale magnetic fields, and second, stronger stratifications 
confine the dominant convective features to a thin outer layer 
where they tend to produce non-axisymmetric field. We had the 
hope that introducing the weakly conducting layer alone would 
promote dipolar dynamo action, allowing zonal winds to de- 
velop in the outer layer while dipolar dynamo action takes place 
in the deeper interior (Stanley and Glatzmaier, 2010; Heimpel 
and Gomez Perez, 201 1). 

The top row of Fig. 6 illustrates how the zonal flow and 
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the poloidal magnetic field change when increasing the density 
stratification in a case with a pronounced weakly conducting 
layer (x, n = 80%) and the bottom row shows the corresponding 
non-magnetic cases. In the latter, increasing density stratifi- 
cation has very significant effect on the retrograde jet, which 
gradually disappears and is mostly gone at N p = 5. This shows 
that increasing N p reduces drastically the amplitude of the zonal 
flows in the conducting part, where the dynamo is generated. 
This asymmetry between prograde and retrograde jets is typi- 
cal for anelastic models with stress-free boundaries and strong 
density stratifications. For N p = 0, the typical three jet structure 
develops with a strong prograde equatorial jet and two flanking 
retrograde jets at higher latitudes (Gastine and Wicht, 2012). 
These flows still have a geostrophic structure and they pene- 
trate the whole shell, in a similar way as the corresponding hy- 
drodynamic solution. The competition between dipolar fields 
and zonal flows leads to the oscillatory solution type (pink area 
in Fig. 2) where either the dipolar field or the zonal flows are 
strong. The remaining prograde jet completely resides within 
the weak conductivity envelope so that the zonal flow/dipole 
competition is not an issue any more. The same is true for 
the confined convective features that additionally promote non- 
axisymmetric magnetic fields at larger stratifications. 

The poloidal magnetic field is now produced by convection 
residing deeper in the shell. This is illustrated by Fig. 7 showing 
the radial profiles of poloidal and toroidal magnetic energies for 
the five solutions depicted in Fig. 6. Up to a stratification of 
N p = 4, poloidal and toroidal energies have similar profiles and 
peak in the outer part of the conducting region. For N p = 5, 
however, the profiles are different with a pronounced focus on 
deeper parts of the shell. For N p < 4 the toroidal field is larger 
than the poloidal, suggesting that the induction mechanism is 
different from the dipole-dominated cases where the poloidal 
field is stronger. 




Figure 7: Radial profile of magnetic energy flux (r 2 E mag ) averaged 
over time. The dashed black line is the location of^ m = 80%. These 
results correspond to the red triangles and red dashed line from Fig. 4. 
The poloidal (dashed lines) and toroidal (dot-dashed lines) compo- 
nents are also shown for N p = 5 and N p = 0, with the corresponding 
colours. The magnetic energy fluxes are normalized by their maxi- 
mum values. 



Gastine et al. (2012) have reported that the multipolar solu- 
tions with stronger zonal flows are dynamos of aQ. or a 2 Q type. 
The Q stands here for the production of toroidal field by global 
zonal wind shear while a stands for poloidal and toroidal field 
production by more local helical structures. It is difficult to dif- 
ferentiate whether the importance of the Q-effect is a reason 
for the multipolar nature of these solutions or arises because 
the large scale Lorentz forces tied to dipolar field contributions 
are missing. 

Following Brown et al. (2011), the Q-effect is given by 

8r\ r I r d0\sme) 
and the mean ohmic diffusion (MD - mean diffusion) by 

MD = ^V-#^(-^l (26) 
r 2 Sin 6» dr\r dr I 

where the overbars correspond to zonal averages. Figure 8 dis- 
plays the axisymmetric azimuthal toroidal field along with the 
Q-effect and the respective mean Ohmic diffusion for the N p = 3 
and the N p -5 solutions already depicted in Figs. 6 and 7. For 
the multipolar solution at N p = 3, the Q-effect plays an impor- 
tant role as demonstrated by the high degree of correlation with 
the azimuthal field over the shell. For the dipole-dominated so- 
lution at N p = 5, however, the Q-effect is only strong in the 
weakly conducting region where it is effectively balanced by 
the mean diffusion. The dynamo mechanism thus changes from 
aQ or a 2 Q to a 2 at stronger stratification. Dynamos of the lat- 
ter type are known to produce dipole-dominated magnetic fields 
(Olson et al., 1999). 

3.3. Ekman number 10~ 5 

The Boussinesq results by Heimpel and Gomez Perez (201 1) 
suggest that the competition between zonal flows and dipo- 
lar field may be less severe at the lower Ekman number 10" 5 
they employed. To explore this, we added eleven simulations 
at E — 10~ 5 and different Rayleigh numbers and stratifica- 
tions, adopting the transition radius Xm - 80% and one case 
of Xm =95%. Note that one Boussinesq case is identical to one 
of the models presented by Heimpel and Gomez Perez (2011) 
(case 50 in Tab. 2). Yellow symbols indicate these runs in both 
Figs. 1 and 2. 

Having values of Roi around 0.01-0.02 and high dipolarity 
they are located on the left side of the dipole-multipole tran- 
sition in Fig. 1. At Roc w 0.035 and N p = 1, there is case 
which has higher value of Rossby number and it is multipolar, 
as expected. Unlike for the larger Ekman number of E — 10~ 4 , 
even Boussinesq or mildly stratified solutions are now dipole- 
dominated and fall into the yellow regime outlined in Fig. 2. 
They are also very stable in time. 

Figure 9 compares two E = 10~ 5 solutions at N p — 1 and 
N p = 3 with the N p — 5 and E - 10~ 4 model already shown in 
Fig. 6. The field is clearly strongly dipolar in all cases, there are 
pronounced equatorial jets and much weaker retrograde jets at 
higher latitudes; both aspects were already reported by Heim- 
pel and Gomez Perez (2011) for their Boussinesq model. At 
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Figure 8: From left to right, azimuthal averages of the toroidal component of the magnetic field, production of the toroidal field by £l-effect and 
ohmic diffusion. The three upper panels correspond to N p = 3, Ra = 2 x 10 7 and the three lower panels to N p = 5, Ra = 4 x 10 7 , both cases belong 
to the red line of Fig. 4. 



both Ekman numbers, the relative amplitude of the retrograde 
jets decreases with N p . Note, however, that the amplitude of the 
dominant equatorial jet is significantly smaller in all lower Ek- 
man number cases. For example, the zonal flow Rossby number 
is Ro zo „ = 9.8 x 1(T 4 in the (E = 1(T 5 , N p = 1) case depicted in 
Fig. 9 but Ro zo „ = 8.0 x 1(T 3 in the (E = 1(T 4 , N p = 5) simula- 
tion. The peak surface velocity of the equatorial jet is around 
Ro eJ = 4.8 x 1(T 2 for Jupiter (Vasavada and Showman, 2005) 
and the E = 10~ 4 simulations come much closer to this value 
than the E = 10~ 5 cases. The (E = 10~ 4 , N p = 5) case illus- 
trated in Fig. 9 reaches indeed Ro e j = 5.3 x 10~ 2 compared to 
Ro ej = 3.7x 10~ 3 for the (£=10- 5 , N p = 1) example. 

The slower jet speeds at the lower Ekman number are not 
caused by Lorentz forces but are already present in the non- 
magnetic simulations. In order to reach more realistic flow am- 
plitudes at E — 10~ 5 the Rayleigh number needs to be increased. 
According to Heimpel and Gomez Perez (2011), however, this 
would once more lead to multipolar magnetic fields. The larger 
jet speed reached in the (E = 10 5 , N p - 3) case shown in Fig. 9 
suggests that increasing the stratification may once more help 
here to retain dipolar dynamo action in the presence of stronger 
equatorial jets. For weak stratification the principle difficulty of 
consolidating strong zonal winds and dipole dominated fields 
remains an issue at the lower Ekman number. 



Elsasser number and magnetic Reynolds number are compa- 
rable in our dipole-dominated E = 10 4 cases of N p > 5 with 
thicker outer low-conductivity layers and the E = 10~ 5 cases 
(see Tab. 1). The Elsasser number is a measure for the ratio of 
Lorentz to Coriolis forces in the Navier-Stokes equation, while 
the magnetic Reynolds number quantifies the ratio of magnetic 
induction to diffusion in the dynamo equation. Similar force 
balance and magnetic balance thus seems to hold at both Ek- 
man numbers with zonal flows playing only a minor role in 
these dominantly a^-dynamos. For a more detailed analyses 
we refer to Heimpel and Gomez Perez (201 1). 

4. Discussion and conclusions 

We have studied dynamo simulations with a conductivity 
profile geared to model the dynamics of the metallic and the 
molecular envelopes of the gas giants in one simulation. The 
goal was to find strong zonal surface jets along with dipole- 
dominated dynamo action which is difficult when the electrical 
conductivity is homogeneous (Gastine et al., 2012). The solu- 
tions then either have strong zonal winds and multipolar mag- 
netic fields or weak zonal winds and dipole-dominated mag- 
netic fields. Sometimes bistability is found where both types of 
solutions coexist at the same parameters. 
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Figure 9: Azimuthal average of the zonal component of the flow. The parameters are: E = 10~ 5 , Xm = 80%, Ra = 4ARa cr for N p = 1 and 
Ra = 3.9 Ra cr for N p = 3. The poloidal fieldlines are plotted on top of the zonal velocity contours, where solid lines are positive and dashed are 
negative values. The dotted lines correspond to r m = 80%. Zonal velocities are given in units of Rossby number calculated by Ro = ti/(f2r ). 



Even thin weakly conducting outer layers already have an 
effect on the dynamo process (Gomez-Perez et al., 2010) and 
reduce the window where dipole-dominated dynamo solutions 
can be found at mild stratifications. The local Rossby number 
that marks the transition between dipolar and multipolar cases 
reduces from Ro c ~ 0.1 to Ro{ * 0.03 tilXm = 95%. The zonal 
flows remain weak in some of these cases and the reason for 
the loss of dipolar dominance remains unclear. Gomez-Perez 
et al. (2010) speculate that the separation of the Ekman and 
Hartmann boundary layers may play a role. 

A sizable weakly conducting outer layer (x m = 80%) al- 
lows stronger zonal winds to develop and the competition with 
the dipolar field becomes even more problematic than in sim- 
ulations with homogeneous conductivity. When the stratifica- 
tion is weak, only multipolar solutions remain where dipole- 
dominated solutions could still be found for homogenous con- 
ductivity. This changes, however, when the stratification is in- 
creased at E = 10~ 4 . At intermediate values, the solutions be- 
come highly time dependent, oscillating between dipolar and 
multipolar configurations. At strong stratifications, highly dipo- 
lar and stable solutions are found. The main reason likely is 
that for the strongly stratified cases explored here the fierce 
equatorial jet is confined to the weakly conducting outer re- 
gion introduced to model the molecular envelope. The flanking 
retrograde jets, however, that cover the higher conducting in- 
ner region where the main dynamo action takes place, remain 
weak. For lower Ekman numbers, dipole-dominated solutions 
are also obtained at lower density stratifications. These solu- 
tions present weaker zonal jets, similarly to previous Boussi- 
nesq studies by Heimpel and Gomez Perez (201 1). 

Compared to Jupiter or Saturn there are too few zonal jets 
in our simulations. Also, the strong decrease in amplitude 
from equatorial to flanking jets necessary to retain dipolar dy- 
namo action is not compatible with the observations, at least for 
Jupiter. Dipolar dynamo action could nevertheless be retained 
for Jupiter if the higher latitude jets were too shallow to become 
important for the dynamo process. 

Another argument against deep reaching winds is that the as- 



sociated strong Q-effect and Ohmic dissipation may also not 
be compatible with Jupiter's observed luminosity (Liu et al., 
2008). First analysis of our results confirm that the Q-effect 
and associated Ohmic dissipation can be significant. However, 
as Glatzmaier (2008) argues, the magnetic field may assume 
a configuration that minimizes the efficiency of the £2-effect. 
Magnetic fieldlines tend to remain perpendicular to the flow 
shear which minimizes Lorentz forces and also the £2-effect. 
Fig. 9 illustrates that the magnetic fieldlines indeed tend to be 
aligned with the rotation axis in the outer part of the shell to re- 
main perpendicular to the strong shear where electrical conduc- 
tivity is still sizable. This helps to minimize Ohmic dissipation 
which nevertheless remains significant in all our simulations 
with strong zonal flows. Further investigation are necessary to 
quantify this effect and extrapolate it to the planetary situation. 

Any problems related to Ohmic dissipation and dipolar dy- 
namo action would be solved if the main zonal jets were con- 
fined to a thin outer envelope where the electrical conduc- 
tivity is negligible (Liu et al., 2008). In our simulations all 
the stronger jets obey the Taylor-Proudman theorem and reach 
through the planet. Shallow jets have been found by Kaspi et al. 
(2009), however, who use a different anelastic approximation 
and a different internal heating mode. Further investigations 
are required to clarify which specific model features influences 
the depth on the zonal jets. 
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Table 2: Summary of the time-averaged results 
* no-slip top boundary 

"a = 25 

""case from Heimpel and Gomez Perez (201 1) with 77 = 0.35 



Model Xm(%) N p E Pirij N r l max Dipolarity SD dip Ro- on Rm Rot A r 

01 80 0.0 11.5 10~ 4 2.0 73 106 4.26xl0" 2 4.73 X 10~ 2 8.73 x 10 -3 115 1.84x 10^ 0.225 3.6 

02 80 0.0 17.2 10~ 4 2.0 73 106 8.13 x 10~ 2 9.46 x 10~ 2 1.41 x 10~ 2 176 3.02 x 10~ 2 0.727 1.4 

03 80 0.0 23.0 10~ 4 2.0 73 106 7.78 x 10~ 2 9.24 x 10~ 2 1.80 x 10~ 2 221 4.13 xl0~ 2 1.238 3.4 
04* 80 0.0 23.0 10~ 4 2.0 73 106 6.94 x 10~ 2 8.43 X 10 -2 5.45 x 10~ 3 180 6.34 xl0~ 2 0.190 1.5 

05 80 0.0 45.9 10~ 4 2.0 81 170 1.53 x 10-' 1.49x10-' 3.09 x 10~ 2 365 7.55 x 10~ 2 2.906 1.0 

06 ~90 00 LL5 lO 4 TO 73 L06 4.89 x 10~ 2 6.07 x 10~ 2 6.86 x 10~ 3 143 2.39 x 10~ 2 0.429 IS 
07* 90 0.0 11.5 10~ 4 2.0 73 106 8.63 x 10-' 2.92xl0~ 2 132xlO- 3 113 2.65 x 10~ 2 1.576 2.9 
08* 90 0.0 23.0 10~ 4 2.0 73 106 1.65 x 10-' 1.50x10-' 4.63 x 10~ 3 237 6.45 x 10~ 2 0.579 2.2 
09d 95 00 LL5 lO 4 To 73 85 8.24 x 10"' 5.75 x 10" 2 1.61 x 10" 3 127 2.72 x 10~ 2 1.836 33" 
09m 95 0.0 11.5 10~ 4 2.0 73 106 1.64x 10-' 9.20xl0~ 2 642xlO~ 3 158 2.66xl0~ 2 0.509 2.1 

10 95 0.0 17.2 10~ 4 2.0 73 106 8.10 x 10~ 2 7.82 X 10~ 2 8.23 X 10~ 3 227 4.46 xlO~ 2 1.295 1.8 

11 95 0.0 23.0 10~ 4 2.0 73 106 7.24 x 10~ 2 7.83 x 10~ 2 944xlO~ 3 282 5.83 x 10~ 2 2.156 1.4 

12 ~80 TO 52 lO 4 To 73 106 1.59 x 10~ 2 3.12 x 10~ 2 8.30 x 10" 3 97 2.24 x 10~ 2 0.230 43 

13 80 1.0 10.3 10~ 4 2.0 73 106 1.40x 10-' 1.40x10-' 2.12 x 10~ 2 220 5.31 x 10~ 2 1.256 2.0 
14d 90 TO 5^2 lO 4 TO 73 85 8.74 x 10" 1 2.83 x 10~ 2 1.89 x 10~ 3 104 3.35 x 10~ 2 1.918 L8 
14m 90 1.0 5.2 10~ 4 2.0 73 85 8.12 x 10~ 2 1.08x10-' 4.72 x 10" 3 125 3.47 x 10~ 2 0.602 1.4 
15* 90 1.0 5.2 10~ 4 2.0 73 85 8.36x 10-' 3.68X10 -2 1.69 x 10" 3 113 3.68x 10^ 2.144 2.0 
16 90 1.0 7.8 10~ 4 2.0 73 106 4.70 xlO" 2 5.80 xlO -2 8.94 xlO -3 208 5.43 X 10" 2 1.526 1.8 
17d 95 LO 41 KF T6 73 85 8.29 x 10~ 2 4.76 x 10~ 2 1.29 x 10 3 91 2.64 x 10" 2 L25l L2 
17m 95 1.0 4.1 10~ 4 2.0 73 85 1.21 x 10-' 9.57xl0~ 2 3.08 xlO -3 106 2.52 xl0~ 2 0.408 1.2 
18d 95 1.0 5.2 10~ 4 2.0 73 106 8.33 x 10-' 4.08 x 10~ 2 1.59 x 10~ 3 123 3.81 x 10~ 2 2.346 2.6 
18m 95 1.0 5.2 10~ 4 2.0 73 106 3.40 x 10~ 2 4.54X10 -2 3.67 x 10~ 3 145 4.16 xlO~ 2 0.770 2.4 

19 95 1.0 7.8 10~ 4 2.0 73 85 3.86 X 10" 2 4.81 x 10 -2 1.00 xl0~ 2 253 6.33 x 10~ 2 1.892 1.1 

~TX) 95 To 2~9 iO 4 T6 73 106 1.94 x 10~ 2 2.08 x 10~ 2 4.64 x 10~ 3 121 3.91 x 10~ 2 0.522 24 

~~ 21 70 3T0 43 lO 4 T6 73 106 1.60 x 10"' 1.77 x 10~' 2.57 x 10~ 2 95 2.26 x 10~ 2 0.204 L6 

22 ~80 3J) 32 UP To 73 85 1.33 x 10~' 9.52 x 10~ 2 1.39 x 10~ 2 78 2.45 x 10~ 2 0.210 31 

23 80 3.0 4.3 10~ 4 2.0 73 106 2.90x 10-' 2.13x10-' 1.98 x 10~ 2 125 4.25 x 10~ 2 0.538 3.9 
24" 80 3.0 4.3 10~ 4 2.0 121 106 2.31 x 10-' 2.07x10"' 2.20x10^ 115 3.76x 10^ 0.448 2.6 
25* 80 3.0 4.3 10~ 4 2.0 73 106 1.67x 10-' 1.42 x 10~' 8.51 x 10 -3 132 5.56 xlO~ 2 0.747 1.6 
26 80 3.0 8.6 10~ 4 2.0 81 170 2.15 x 10-' 1.84x10-' 4.04 x 10~ 2 276 8.57 x 10~ 2 2.372 1.2 

~T7 90 3J) 32 lO 4 To 73 85 5.22 x 10" 2 7.13 x lO 2 9.37 x 10~ 3 133 4.16 x 10~ 2 0.633 24 

"38 95 To 32 lO 4 To 73 106 3.31 x 10 -2 6.71 x 10 -2 7.57 x 10~ 3 165 6.00 x 10~ 2 0.956 JI 

~~29 70 40 8^8 HP To 81 170 1.64 x 10" 1 1.37 x 10" 1 4.63 x 10" 2 137 3.01 x 10" 2 0.422 L7 

~T6^ 80 40 44 TO 81 170 1.21 x 10~' 1.22 x 10 _1 1.52 x 10 2 67 1.71 x 10~ 2 0.092 33 

31 80 4.0 5.5 lO" 4 2.0 81 170 5.04x 10-' 3.05x10-' 2.01 x 10~ 2 98 3.27 x 10~ 2 0.381 2.3 

32 80 4.0 8.8 10~ 4 2.0 81 170 2.36x 10-' 1.85 x 10 _1 3.04x10^ 176 5.45 X 10~ 2 1.142 1.2 
~T3 95 40^ 53 W 4 TO 81 170 6.73 x lO" 4 8.86 x 10" 4 9.73 x 10" 3 198 8.57 x 10" 2 1.626 L3 
~T4^ 80 10 74 W 1 TO 81 170 8.97 x 10~' 2.12 x 10^ 2 1.44 x 10 2 71 3.33 x 10^ 2 0.829 2~9 

35" 80 5.0 7.4 10~ 4 2.0 97 170 5.29 x 10-' 2.44x10-' 2.15xl0~ 2 69 2.11 xl0- 2 0.240 1.6 

36d 80 5.0 9.3 10~ 4 2.0 81 170 8.73 x 10-' 1.64xl0~ 2 1.58 xlO -2 94 4.60x 10^ 1.082 1.3 

36m 80 5.0 9.3 10~ 4 2.0 97 170 2.93 x 10-' 2.28 x 10 _1 2.29 xlO" 2 104 4.34x 10^ 0.616 1.6 

37 90 5l) 74 1(P TO 97 170 3.15 x 10~ 2 ToTxW 1 1.02 x 10" 2 145 5.39 x 10^ 2 1.135 LT 

38 90 5.0 9.3 10~ 4 2.0 97 170 3.47 x 10~ 2 3.29 xlO -2 9.75 X 10 -3 180 7.43 X 10" 2 1.689 0.7 
"T9 80 53 91 lO 4 T6 97 170 9.14 x 10 _1 1.53 x 10~ 2 1.39 x 10^ 2 76 3.44 x 10~ 2 0.820 L4 

~~ 40 80 00 12 W 5 LO 81 170 8.56 x 10" 1 2.40 x 10" 2 5.01 x lO" 4 79 7.44 x 10" 3 0.319 LO 

41 80 0.0 12.1 10~ 5 1.0 81 170 8.45 x 10-' 2.12X10 -2 1.01 x 10~ 3 122 1.37x 10^ 1.090 0.9 

95 00 5T8 10~ 5 LO 81 133 8.51 x 10" 1 1.16 x 10" 2 3.21 x lO" 4 97 7.05 x 10~ 3 0.209 IT 

~^3 80 LO 44 lO 3 LO 81 170 9.25 x 10~' 8.82 x lO" 3 1.55 x 10 3 96 1.32 x 10" 2 0.538 TT 

44 80 1.0 6.7 10~ 5 1.0 81 170 9.27 x 10-' 1.21 x 10~ 2 3.67 x 10" 3 162 1.99 x 10~ 2 1.071 1.1 

45 80 1.0 11.1 10~ 5 1.0 97 170 6.13xl0" 2 8.68x 10^ 7.19xl0~ 3 316 3.62 X 10" 2 0.763 0.2 
"36 80 TO T6 lb~ 5 LO 81 i70 9.22 x IO-' 7.91 x 10~ 3 1.90 x lO" 3 74 1.31 x 10~ 2 0.432 08 

47 80 2.0 5.4 10~ 5 1.0 81 170 8.84 x 10~' 3.82x10^ 4.79x10^ 149 1.82x 10^ 0.485 1.5 

Continued on next page 
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Model 






Ra 

Rd cr 


E 


Pvrii 




^max 


Dipolarity 


SD dip 


Ro zo „ 


Rm 


Ro t 


A 


T 


48 
49 


80 
80 


3.0 
3.0 


2.9 
3.9 


io- 5 

IO" 5 


1.0 
1.0 


81 
81 


170 
170 


9.50 x 10-' 
8.74 x IO" 1 


1.16 x 10~ 2 
3.05 x IO" 2 


4.58 x IO" 3 
4.57 x IO" 3 


100 
144 


1.43 x IO" 2 
2.42 x IO" 2 


0.348 
0.853 


1.5 
1.0 


50"* 


80 


0.0 


7.2 


io- 5 


3.0 


121 


170 


8.05 x 10-' 


1.46 x 10~ 2 


1.74 x IO" 3 


306 


1.20 x IO- 2 


2.279 


0.2 
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